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SHOCK  WAVE  PROPAGATION  IN  AN  INHOMOGENEOUS  MEDIUM 
USING  FINITE  DIFFERENCES 

Section  1 

INTRODUCTION 

* 1-4 

Finite-difference  fluid  equations  have  been  known  to  pro- 
vide valid  solutions  to  problems  containing  shocks  when  the  physically 
correct  conservation  variables  in  conservation  form  are  used.  In  this 
study,  we  examine  the  problem  of  shock  propagation  in  an  inhomogeneous 
medium  with  exponentially  varying  density.  The  one-dimensional  self- 
similar analytic  solution  to  this  problem  will  be  compared  to  various 
numerical  solutions  using  different  algorithms.  This  problem  consti- 
tutes a more  severe  numerical  test  for  shock  propagation  than  the  shock 
wave  in  an  homogeneous  medium  and  therefore  serves  as  a better  test  of 
numerical  algorithms. 

In  particular  we  find  that  care  must  be  taken  in  the  use  of  forms 
of  the  hydrodynamic  equations  which  do  not  express  physical  conservation. 
For  non-conservation  formulations  of  the  energy  equation  an  artificial 
viscosity  must  be  introduced,  not  only  to  provide  the  necessary  stability, 
but  also  to  provide  shock  heating.  The  magnitude  of  firs  artificial  vis- 
cosity to  obtain  best  shock  results  depends  on  the  grid  size  and  the 
problem  type.  There  is  no  simple  way  to  obtain  this  optimal  viscosity 
for  problems  where  the  solution  is  not  known  in  advance. 
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In  addition,  we  find  that  Flux- Corrected  Transport  (FCT)  has 
several  properties  which  make  it  more  flexible  and  effective  for 
shock  calculations.  While  the  comparison  between  different  energy 
formulations  has  been  made  easier  through  the  use  of  the  FCT  scheme, 
the  results  hold  for  any  finite  difference  algorithm  and  in  particular 
they  will  be  shown  to  hold  using  the  Lax-Wendroff  scheme  as  well. 


Note:  Manuscript  submitted  August  22,  1975. 
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T he  problem  and  its  analytical  solution  are  well  described 

g 

by  Zeldovich  and  Raizer.  We  outline  here  their  self-similar  solution 

as  well  as  the  conditions  under  which  it  holds.  The  solution  is  given 

both  for  the  increasing  and  decreasing  density  cases.  The  medium  is 

defined  by  its  characteristic  length  - or  scale  height  - which  is  equal 

to  the  e-folding  distance.  The  analytic  solution  is  given  as  a function 

of  the  similarity  variable  £ which  is  equal  to  — ^ — where  X is  the 

shock  location  and  x is  the  Eulerian  coordinate  location.  The  velocity 

of  the  shock  front  is  given  by  D where  a - a coefficient  which 

depends  only  on  the  specific  heat  ratio  V - is  determined  by  the  solution 
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of  a boundary-value  problem.  For  the  increasing  density  problem 
with  7 = 2 (a  = 1. 5),  the  self  similar  analytic  solution  behind  the  shock  is: 

£4  p0o;)  (1  + 2 f)‘5/2 
yfl  “2  ^ pc(x)  a + 2 f)"3/2 

yh.  « f « - « 

density  immediately  in  front  of  the  shock, 
characterized  by 
P0(X)  e‘f 

0 

0 

The  initial  conditions  for  our  calculations  were  also  taken  after  the 
self -similar  solution  had  developed.  The  only  restriction  on  the 
application  of  the  analytic  solution  is  that  the  initial  pressure  and 
temperature  be  taken  equal  to  0.  While  this  condition  is  difficult  to 
satisfy  rigorously  in  a numerical  calculation,  the  temperature  and 


P = 
P = 

u = 

where  PQ(X)  is  the  ambient 
The  unperturbed  medium  is 

P = 

P = 
u = 


2 


the  pressure  of  the  ambient  gas  have  been  chosen  to  be  very  small  so 

that  the  pressure  and  temperature  ratios  across  the  shock  are  very 
4 

large  (5x  10  typically  for  the  pressure  ratio).  This  corresponds  to 
the  infinite  Mach  number,  strong  shock  limit  required  for  the  self- 
similar solution. 

Three  different  variables  have  been  used  for  the  energy 
equation,  while  the  continuity  and  momentum  equations  have  been 
treated  in  their  usual  conservative  form.  The  numerical  algorithms 
used  treat  all  conservative  terms  in  conservative  finite-difference 
form.  Use  has  been  made  of  a Lax-Wendroff  scheme  (LW)  and  a 
Flux-Corrected  Transport  (FCT)  scheme.  The  Lax-Wendroff  scheme 
uses  the  two-step  Richtmeyer  form  coupled  with  a Von  Neumann  vis- 
cosity to  provide  the  additional  stability  and  viscous  heating  needed 
near  the  shocks.  The  FCT  scheme  makes  use  of  artificial  viscosity 
(Von  Neumann  type)  when  the  energy  equation  is  cast  in  terms  of  the 
pressure  or  temperature  variable,  to  provide  the  viscous  heating 
otherwise  lacking  in  the  shocks.  Artificial  viscosity  is  not  needed 
in  the  strict  conservative  formulation. 

The  first  set  of  equations  used  solves  for  the  mass,  momen- 
tum, and  total  energy  per  unit  volume  and  is  in  strict  conservative 
form.  The  second  set  is  based  on  mass,  momentum,  and  pressure 
whereas  the  last  one  uses  mass,  momentum,  and  temperature.  The 
total  energy,  pressure,  and  temperature  equations  are  respectively: 

- v • (p  + q)  v (1) 

- (?-  l)(p  + q)  V • v (2) 

- [(r-  2)  T + (y-i)|]  v • v 


dE 

v-vE  = 


ap 

9F+2'I  p * 


<9T  T _ 

77  * 
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(3) 


(4) 


q = .Pb(6x2>|-!ii. 

12  n 

where  E = p(c+  ~ v ) and  e = for  a perfect  gas. 

The  conservation  form  of  the  difference  equations  is  a necessary 
but  not  a sufficient  condition  to  guarantee  the  correct  weak  (discontin- 
uous) solution  to  the  ideal  hydrodynamic  equations.  It  is  possible  to 
get  different  weak  solutions  by  using  different  equivalent  forms  of  the 
partial  differential  equations.  The  jump  conditions  depend  on  the 
form  of  the  equation  and  we  must  use  physical  reasoning  to  determine 
if  these  jump  conditions  make  sense » In  dealing  with  physical  laws 
we  usually  try  to  write  the  equations  in  a conservation  form  which 
implies  actual  physical  conservation.  In  hydrodynamics  only  the 
total  energy  equation  combined  with  the  mass  and  momentum  equation 
provides  a consistent  set  of  equations  for  physically  conserved  quantities. 

This  problem  of  correctly  calculating  the  weak  solutions  can  be 
avoided  by  introducing  viscosity  and  seeking  genuine  (continuous) 
solutions.  With  viscosity  the  hydrodynamic  equations  do  not  have  truly 
discontinuous  solutions  and  all  forms  of  the  energy  equation  should  be 
equivalent.  Artificial  viscosities  - and  not  real  viscosity  - have  to  be 
added  to  provide  the  necessary  heating  because  real  viscosity  acts 
over  a few  mean  free  paths  which  are  usually  much  smaller  than  the 
grid  size.  The  effect  of  a real  physical  viscosity  will  always  be  too 
small  unless  the  mean  free  path  and  cell  size  are  comparable. 

However,  the  difficulty  with  artificial  viscosity  is  that  the  shock 
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profile  must  be  spread  over  several  cells  for  stability  regardless 
of  the  actual  physical  thickness  of  the  shock  and  thus  inter- 
actions of  the  shock  with  inhomogeneous  media  will  not  be  accurate 
unless  the  background  varies  slowly.  Further,  since  there 
is  little  physical  basis  for  the  artificial  viscosity,  we  can  only  hope 
that  it  will  produce  the  correct  heating  in  non-conservative  formulations. 

In  this  regard  the  FCT  finite  difference  formulation  seems  to 
have  a distinct  advantage  in  that  its  shock  profile  spreads  over  only 
a couple  of  cell  widths  almost  independent  of  the  coefficient  used  for 
the  artificial  viscosity.  The  velocity  gradient  used  in  the  Von  Neumann 
viscosity  is  nearly  independent  of  the  coefficient  and  thus  any  amount 
of  heating  that  is  desired  can  be  achieved  by  raising  the  viscosity 
coefficient.  In  most  other  schemes  the  larger  coefficient  tends  to 
smooth  out  the  shock  profile  which  reduces  the  velocity  gradient 
and  somewhat  offsets  the  effect  of  increased  coefficient. 

It  is  the  failure  to  compute  accurately  the  dissipation  mech- 
anism which  converts  kinetic  energy  to  thermal  energy  that  leads  to 
a failure  to  conserve  energy  in  the  temperature  and  pressure  formu- 
lations and  hence  gives  incorrect  results  for  the  shock  dynamics. 

Since  the  total  energy  equation  is  in  divergence  form  whether 
the  viscosity  terms  are  included  or  not,  conservation  of  energy  is 
automatically  guaranteed  when  a conservative  difference  scheme  is 
applied  to  it.  In  the  remainder  of  the  paper  we  will  show  the  results 
of  several  test  calculations  demonstrating  these  ideas.  Section  Two 
shows  the  results  for  a shock  propagating  into  an  exponentially 
increasing  density  medium.  The  results  of  the  different  formulations 
are  compared  for  several  values  of  the  artificial  viscosity  parameters 
and  grid  sizes.  In  Section  Three  the  results  for  the  decreasing  density 
case  are  shown  and  in  Section  Four  the  conclusions  that  can  be 
drawn  from  this  study  are  made. 
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Section  2 


RESULTS  FOR  THE  INCREASING  DENSITY  CASE 


For  a y=2  gas  the  self- similar  solution  is  completely  analytic. 
From  the  Rankine-Hugoniot  relations,  we  expect  the  density  jump 
across  the  shock  to  be  equal  to  3.  Figures  1 and  2 show  density  profiles 
for  the  energy  and  temperature  equations  respectively  for  the  LW 
scheme  after  a time  t = 350  St.  At  that  time,  the  shock  has  moved  over 
a distance  equal  to  1. 3 A.  Three  different  values  of  the  viscosity 
coefficient  have  been  used  and  the  effect  of  non-conservation  is  shown 
clearly  from  these  two  figures.  In  the  total  energy  formulation,  the 
value  of  b affects  mainly  the  stability  of  the  solution  (and  the  amplitude 
of  the  ripples  behind  the  shock);  in  the  temperature  formulation,  it 
changes  the  speed  of  propagation  of  the  shock  significantly.  For  the 
latter  equation,  the  larger  the  viscosity  coefficient,  the  more  the 
viscous  heating  and  the  better  the  agreement  between  the  numerical 
solution  and  the  analytic  solution.  Note,  however,  that  the  peak 
density  behind  the  shock  decreases  with  the  artificial  viscosity 
coefficient  b and  thus  the  density  profiles  cannot  be  taken  as  the 
only  criteria  of  good  numerical  solu^ons. 

The  integrated  total  and  thermal  energies  are  shown  in 
Figure  3 for  some  of  these  cases.  This  figure  shows  that  energy 
conservation  is  better  achieved  by  using  the  total  energy  formulation. 
Note  that  while  for  the  correct  energy  conserving  formulation  the 
integrated  thermal  energy  is  increasing  with  time,  in  the  case  of 
the  temperature  formulation  it  is  actually  decreasing.  As 
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we  have  seen,  this  loss  of  thermal  energy  reduces  the  driving 
force  of  the  shock  and  results  in  the  shock  lagging  the  correct 
solution.  When  the  temperature  equation  is  used,  energy 
conservation  is  improved  by  increasing  the  magnitude  of  the  artificial 
viscosity  coefficient.  A limiting  case  occurs  when  the  b coefficient 
in  q is  sufficiently  large  to  make  the  artificial  pressure  greater  than 
the  real  pressure  in  the  shock  front.  In  that  case,  too  much  heat  is 
generated  in  the  shock  and,  for  example  by  using  b=4  in  the  FCT 
temperature  scheme  the  shock  moves  faster  than  its  analytic 
counterpart  by  1%.  The  total  integrated  energy  also  increases  above 
its  correct  value  for  this  case  and  the  difference  reaches  1. 8%  at  the 
end  of  the  run.  Note  also  that  for  these  cases  the  stability  requires 
a smaller  time  step  since  both  the  temperature  and  the  shock  velocity 
take  larger  values  than  in  a real  shock. 

Another  limiting  case  is  shown  in  Figure  4 where  q has  been 
set  equal  to  0 for  both  the  temperature  and  ;ne  pressure  equation 
formulation.  The  algorithm  used  in  this  case  is  FCT  since  LW  would 
be  unstable.  The  shock  lags  behind  its  exact  solution  to  a much  larger 
extent  than  shown  previously  and  36%  of  the  total  energy  is  lost  using 
this  formulation. 

It  was  investigated  whether  these  results  depend  on  the  form 
of  the  artificial  viscosity  coefficient.  In  Figure  5,  use  has  been  made 
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of  the  Tyler  viscosity  coefficient  with  the  LW  scheme  using  the 
temperature  formulation.  The  density  profiles  look  smoother  than 
for  the  Von  Neumann  viscosity  but  the  energy  conservation  is  not  as 
good  and  consequently  the  shock  location  is  also  worse.  Once  again, 
smooth  density  profiles  do  not  constitute  a complete  criterion  for 
good  numerical  solutions.  The  shock  dynamics,  which  can  be  checked 
only  against  an  exact  solution,  must  also  be  considered. 
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Variations  of  the  Results  with  Grid  Size 


Results  shown  previously  have  been  obtained  with  a fixed 
grid  size  corresponding  to  a spatial  resolution  of  40  grid  points  per 
scale  height.  In  practice,  the  resolution  is  often  much  coarser  so 
the  influence  of  the  spatial  resolution  on  the  results  is  now  investi- 
gated. Figure  6 shows  the  results  for  the  density  profiles  when  the 
grid  size  3x  is  multiplied  by  4 so  the  scale  height  is  made  up  of  10 
grid  points;  the  density  profile  in  the  shock  broadens  and  the  peak 
value  of  the  density  just  behind  the  shock  decreases.  For  this 
kind  of  spatial  resolution,  all  equations  have  difficulty  in  simulating 
the  presence  of  a strong  shock  and  in  fact  look  similar.  The  two-cell- 

Q 

wide  flat  top  on  the  density  profile  is  characteristic  of  FCT  . The  energy 
equation,  although  showing  a reduced  density  ratio  across  the  shock, 
still  approximates  fairly  well  the  shock  location  and  yields  energy 
conservation.  The  two  other  energy  equation  formulations  gain 
energy  by  10%  to  16%  with  the  pressure  equation  more  nearly  conserving 
energy.  These  results  contrast  with  the  finer  resolution  cases  for 
which  the  same  values  of  b lead  to  a loss  of  energy.  Although  the  effect 
of  grid  size  is  supposed  to  be  scaled  out  of  the  problem  by  the  form  of 
the  artificial  viscosity  used,  in  effect  when  the  non-conservative 
formulations  are  used  energy  conservation  and  shock  location  are 
altered  by  changes  in  grid  resolution  even  when  the  same  coefficient 
of  artificial  viscosity  is  used.  Clearly,  for  this  resolution, 
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information  about  the  shock  has  been  mostly  lost  for  both  temper- 
ature and. pressure  equations  and  suggests  that  even  10  grid  points 
per  scale  height  with  the  energy  equation  represents  a minimum 
resolution  in  order  to  provide  a meaningful  solution. 

Figure  7 shows  density  profiles  for  the  temperature  equation 
when  the  spatial  resolution  has  been  increased  to  80  cells  per  scale 
height.  The  shock  width  decreases  while  the  density  jump  ratio 


increases  for  the  same  value  of  the  parameter  b in  the  artificial 
viscosity  coefficient. 

From  the  results  presented  for  the  increasing  density  case, 

it  is  apparent  that  only  the  total  energy  equation  formulation  yields  a 

correct  result  in  w^ich  the  shock  location  does  not  depend  strongly 

on  the  viscosity  coefficient  or  on  changes  in  the  grid  size.  For  the 

other  energy  equation  formulations,  although  it  is  possible  to  find  an 

optimum  value  for  the  viscosity  coefficient  in  each  specific  case, 

this  value  is  not  independent  of  changes  in  grid  size  or  problem 

parameters.  The  total  energy  equation  is  thus  superior  in  all 

respects  to  the  non- conservative  forms.  Although  this  result  is 
4 

already  known  , attempts  to  use  the  non-conservative  formulation 
with  (or  even  without)  artificial  viscosity  for  viscous  heating  in  shocks 
have  been  madQ  repeatedly.  Further,  this  work  has  allowed  us  to 
quantify  this  notion  for  a specific  case  by  estimating  the  error  made 
when  a non- conservative  form  is  used. 
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Section  3 


RESULTS  FOR  DECREASING  DENSITY  CASE 


In  case  of  an  exponentially  decreasing  density  medium, 
the  ’’analytic”  self-similar  solution  does  not  exist  and  has  to  be 
replaced  by  the  numerical  solution  of  the  following  ordinary  differential 
equations: 

1 . y-hl  - 2-ghT1/y --(2 “<x)/ay 
dp  _ y+1  y+l  y ? ^ 

4 = 'T'~ Hr(l7y 

1 TF  P 71 


p ry  7?(2/«)+7-l 


= 1 


— du 

n + cctj  -g- 


1 

dij 


« 

dTj 


7+1  ~ AC 


where  u,  p,  and  p is  the  solution  behind  the  shock.  Since  the  solution 
to  these  equations  is  well-behaved,  monotoniq  and  does  not  involve 
computing  the  shock  wave,  good  precision  can  be  obtained  and  a 
meaningful  comparison  can  be  performed  with  the  results  of  the 
partial  differential  equations  used  previously  for  the  "numerical” 
solution.  The  specific  heat  ratio  y was  chosen  to  be  equal  to  7/5 
for  this  case  (a  = 5. 45)  and  results  are  summarized  briefly  below. 

Figures  8 and  9 show  the  density  profiles  for  a very  strong 
shock  propagating  in  a decreasing  density  medium  for  the  FCT  and 
LW  schemes  respectively.  At  the  time  it  is  shown  (t  = 200 St)  the 
shock  has  traveled  approximately  a distance  equal  to  1. 2 A.  The  grid 
size  is  the  same  as  that  used  in  Figures  1-5.  Note  that  this  time  the 
shock  is  accelerating. 
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Several  interesting  features  may  be  noted  from  these  graphs. 
First,  using  an  artificial  viscosity  coefficient  h= 2 results  in  the 
shock  propagating  too  fast  for  both  the  FCT  and  LW  schemes.  This 
is  in  contrast  to  the  results  of  the  increasing  density  case  shown  in 
Figures  2,  5,  and  7 where  a coefficient  b=2  resulted  in  too  slow  a 
shock  propagation  A.gain  for  the  decreasing  density  case  a coefficient 
b=0. 5 results  in  the  shock  lagging  behind  the  analytic  solution.  The 
b=l  coefficient  for  the  LW  method  gives  reasonable  shock  location 
but  shows  significant  ripples  behind  the  shock.  For  comparison  the 
solution  obtained  with  the  conservation  energy  formulation  is  shown 
in  Figure  8 and  again  this  shock  location  result  is  in  agreement  with 
the  analytic  solution. 

Integrated  total  and  thermal  energies  can  be  computed  as 
previously.  The  differences  are  much  smaller  for  this  decreasing 
density  case  (of  the  order  of  3%)and  this  can  easily  be  explained  by 
the  fact  that  since  the  shock  propagates  in  a region  of  decreasing 
density,  the  energy  carried  by  the  shock  represents  a decreasingly 
smaller  fraction  of  the  total  initial  energy.  Thus,  while  the  energy 
is  conserved  to  a much  better  accuracy  than  previously  the  error  in 
the  shock  dynamics  is  not  reflected  as  much  in  the  total  energy  con- 
servation and  energy  conservation  is  a less  useful  check  on  accuracy. 

Another  test  shows  the  results  obtained  with  a pressure 
equation.  In  general,  the  use  of  the  pressure  equation  leads  to  the 
same  kind  of  results  as  those  obtained  from  the  temperature  equation. 
More  specifically,  because  it  involves  the  conservation  of  thermal 
euergy,  the  results  it  yields  lie  between  those  given  by  the  temperature 
and  the  total  energy  equation.  The  artificial  viscosity  coefficient  is 
similar  to  that  used  previously  and  we  see  in  this  particular  example 


ll 


that  the  pressure  equation,  although  predicting  quite  accurately  the 

shock  location  using  FCT,  oscillates  much  more  than  the  temperature 

equation.  Oscillations  in  Figs.  8 and  9 may  have  different  origins. 

They  are  clearly  characteristic  oscillations  behind  a shock  for  the 

LW  algorithm.  As  for  the  FCT  scheme,  the  fluctuations  also 

0 

seem  to  be  characteristic  , but  they  are  somewhat  weaker  terraces 
in  which  only  the  derivative  oscillates. 
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Section  4 


CONCLUSIONS 


In  this  study,  it  has  been  shown  that  the  numerical  results 
obtained  for  shock  speeds  and  instantaneous  profile  in  an  exponentially 
varying  density  medium  can  differ  largely  due  to  the  choice  of  energy 
equation  and  spatial  resolution.  By  comparison  with  an  analytic 
solution,  it  has  been  shown  that  only  the  conservative  energy  equation 
is  reliable.  Even  in  this  best  case,  a fairly  fine  spatial  resolution  is 
needed  in  order  to  derive  accurate  results. 

The  inclusion  of  some  artificial  viscosity  is  necessary  not 
only  for  stability  but  to  produce  the  necessary  shock  heating  in  the 
case  of  the  temperature  and  pressure  formulations.  By  suitable 
adjustment  of  the  coefficient  of  artificial  viscosity  one  can  obtain  a 
wide  range  of  shock  profiles  and  shock  heating  and  achieve 
near  conservation  and  therefore  good  solutions.  However, 
it  was  found  that  there  is  no  unique  way  to  choose  this  coefficient 
and  the  precise  value  to  achieve  conservation  depends  both  on  the 
grid  size  and  the  nature  of  the  problem. 

The  FCT  algorithm  does  not  require  artificial  viscosity  for 
stability  and  maintains  a steep  profile  rather  independent  of  the  value 
of  artificial  viscosity.  Thus,  if  the  temperature  or  pressure  equation 
must  be  used,  FCT  gives  more  flexibility  in  achieving  the  correct 
amount  of  heating  in  the  shock  front.  In  addition,  in  the  case  of  the 
total  energy  formu’ation,  the  FCT  scheme  requires  no  artificial 
viscosity  at  all,  removing  an  additional  restriction  on  the  time  step 
and  allowing  larger  time  steps  to  be  used. 


13 


gjajgg. 


ACKNOWLED  GMENTS 

We  are  indebted  to  Dr.  S.  Zalesak  for  lending  us  an  updated 
version  of  FCT.  We  would  like  to  acknowledge  many  useful 
discussions  with  Drs.  D.  Book  and  N.  Winsor.  Finally,  we  thank 
Dr.  J.  Boris  for  his  critical  remarks,  helpful  suggestions,  and 
his  careful  reading  of  the  manuscript. 

This  work  was  supported  by  the  Defense  Nuclear  Agency. 


REFERENCES 


1.  P.  D.  Lax,  Comm.  Pure  Appl.  Math.  7,  159  (1954). 

2.  P.  D.  Lax,  Comm.  Pure  Appl.  Math.  1£,  537  (1957). 

3.  P.  D.  Lax  and  B.  Wendroff,  Comm.  Pure  Appl.  Math.  13, 

217  (1960).  ~ 

4.  J.  Gary,  Math,  of  Comp.  18,  1 (1966). 

5.  J.  P.  Boris  and  D.  L.  Book,  J.  Comput.  Phys.  11,  38  (1973). 

6.  Y.  B.  Zeldovich  and  Y.  P.  Raizer,  Physics  of  Shock  Waves 

and  High  Temperature  Phenomena,  Academic  Press,  Vol  n, 
p.  852-863  (1S56J! 

7.  L.  D.  Tyler  and  M.  A.  Ellis,  "The  Tshok  code:  Lax  Version, M 
SC-TM-70-153,  Sandia  Laboratories  (1970). 

8.  D.  L.  Book,  J.  P.  Boris  and  K.  Hain,  Generalization  of  the 
Flux  Corrected  Transport  Technique,  NRL  Memorandum 
Report  3021  (1975). 


TOTAL  ENERGY  EQUATION- LW  ALGORITHM  - 


X (ARBITRARY  UNITS) 

Fig.  1 — Shock  density  profiles  for  shock  propagating  in  the  increasing  density 
direction.  Total  energy  equation  formulation  with  Lax-Wendroff  (LW)  algorithm. 
The  shock  was  located  at  x = 0 at  t = 0. 
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TEMPERATURE  EQUATION- LW  ALGORITHM  - 


Fig.  2 — Shock  density  profiles  corresponding  to  the  case  of 
Fig.  1 for  temperature  equation  formulation 
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Fig.  3 — Spatially  integrated  total  and  thermal  energies  as  a 
function  of  time  for  cases  of  Figs.  1 and  2 


VARIOUS  ENERGY  EQUATIONS 
-b=0  -FCT  ALGORITHM- 


X (ARBITRARY  UNITS) 

Fig.  4 — Comparison  of  the  different  energy  equation  formulations  (energy  E, 
temperature  T,  pressure  P)  without  any  artificial  viscosity  using  Flux-Corrected 
Transport  (FCT) 


18 


TEMPERATURE  EQUATION -LW  ALGORITHM  - 


X (ARBITRARY  UNITS) 

Fig.  5 — Same  as  Fig.  2 using  Tyler’s  form  of  artificial  viscosity 
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f.  (ARBITRARY  UNITS) 


Kg.  6 — Influence  of  grid  size.  Grid  size  is  four  times  as  large  as  for  previous 
figures.  Shock  density  profiles  are  shown  for  various  en  rgy  equations  formu- 
lations using  FCT  and  for  the  total  energy  equation  using  LW(Elw). 
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TEMPERATURE  EQUATION -LW  ALGORITHM - 


X (ARBITRARY  UNITS) 

Fig.  7 — Influence  of  grid  size.  Grid  size  is  half  that  used  in  Figs. 
1 through  5.  Temperature  equation  formulation  using  LW. 
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TEMPERATURE  AND  TOTAL  ENERGY  EQUATION 
- FCT  ALGORITHM - 


L* 


X (ARBITRARY  UNITS) 

Fig.  8 — Sb^ck  density  profiles  for  shock  propagating  in  the  decreasing  density 
region.  Temperature  and  energy  equation  formulations  using  FCT.  The  shock 
was  located  at  x = 0 at  t = 0. 
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